The Response of Metal Rich Gas to X-Ray Irradiation from a 
Massive Black Hole at High Redshift: Proof of Concept 
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ABSTRACT 

Observational studies show that there is a strong hnk between the formation and evolution 
of galaxies and the growth of supermassive black holes (SMBH) at their centers. However, the 
underlying physics of this observed relation is poorly understood. In order to study the effects of 
X-ray radiation on the surroundings of the black hole, we implement X-ray Dominated Region 
(XDR) physics into Enzo and use the radiation transport module Moray to calculate the radiative 
transfer for a polychromatic spectrum. In this work, we investigate the effects of X-ray irradiation, 
produced by a central massive black hole (MBH) with a mass of M = 5 x 10"* Mq, on ambient 
gas with solar and zero metallicity. We find that in the solar metallicity case, due to high opacity 
of the metals, the energy deposition rate in the central region (< 20 pc) is high and hence the 
temperatures in the center are on the order of 10^~^ K. Moreover, due to the cooling ability and 
high intrinsic opacity of solar metallicity gas, column densities of 10^'' cm~^ are reached at a 
radius of 20 pc from the MBH. These column densities are about 3 orders of magnitudes higher 
than in the zero metallicity case. Furthermore, in the zero metallicity case an X-ray induced H 
II region is formed already after 5.8 Myr. This causes a significant outflow of gas (~ 8 x 10^ M©) 
from the central region, with the gas reaching outflow velocities up to ^ 100 km s~^. At later 
times, ~ 23 Myr after we insert the MBH, we find that the solar metallicity case also develops 
an X-ray induced H ii region, but delayed by ^ 17 Myr. 

Subject headings: galaxies: active galaxies; formation galaxies: nuclei 



INTRODUCTION 



Most galaxies today are thought to host super- 



massive black holes (SMBHs) in their centers (Ko- 
mendy fc Richstone||1995| [Richstone et al.lllOQS 



Kormendy & Gebhardt 2001). Dynamical stud 



Magorrian et al. (1998)) or the velocity dispersion 
(Mbh oc ctP, p ~ 4, |Ferrarese & Merritt| (|2000|); 



Gebhardt et al. (2000)) of the stellar component 



of the host galaxy. More recent studies have shown 
that the value of p varies between 4-6 depending 



ies of central black holes in nearby galaxies reveal 
that there is a relation between the central black 
hole mass and the bulge mass (Mbh /Mb ^ 10"'^, 



on the type of the galaxy ( Graham et al. 2011 Mc 
Connell et al.|[20TT). These relations suggest that 
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there is a strong link between the formation and 
evolution of galaxies, and the growth of SMBHs 
at their centers. 

The black hole masses of quasars observed at 
high redshifts {z > 6) are on the order of 10^ M© 
(iFan et al.||2003l |2006| iKurk et al.||2007|). These 



SMBHs are thought to form (in less than 1 bil- 



lion years) through the accretion of gas onto seed 
black holes of masses between 10^ - 10^ Mp). The 



elling ([Kauffmann fc Haehnelt|[2000{ [Wyithe &: 



origin of these seed black holes (Bromm & Loeb 



2003 Spaans & Silk 2006) is still an open ques- 



tion. However, the existence of quasars with in- 
ferred black hole masses in excess of 10^ Mq im- 
plies that black hole growth should be dominated 
by high accretion rates (> 0.1 Eddington) during 
their quasi-stellar object (QSO) phase. 

The amount of gas available for the SMBHs 
to accrete depends on the mergers or interactions 
with galaxies, which drive large quantities of gas 



into galaxy centers (Barnes & Hernquist 19911, 



star formation and feedback processes. Obser- 
vatio nal (Merloni 20041 and theoretical (Hopkins] 
et al. 2006) studies indicate that there is an evo- 



lution in the relation between the black hole mass 
and the host stellar mass such that at high red- 
shifts (z > 6), the ratio between the black hole 
mass and the host stellar mass is on the order of 
0.02-0.1, which is more than one order of magni- 
tude larger than the ratio measured in local galax- 
ies ( [Peng et al.||2006[ [Murray et al.||2009l ). 

The growth of SMBHs during a quasar phase 
can be traced by the luminosity function of QSOs 
as a function of redshift ( Soltan|[l982 ). The peak 
activity of luminous QSOs occurs at z ~ 2 — 3, 
where the majority of the most massive galax- 
ies were also forming most of their stars. On the 
other hand, the peak activity of lower luminosity 
active galactic nuclei (AGN) occurs at z ~ 0.7— 1. 
X-ray selected AGN serve as a population to ro- 
bustly determine the luminosity function and evo- 
lution of unobscured AGN. However, the obscured 
fraction of AGN depends on X-ray luminosity as 
well as redshift. Therefore, this obscured frac- 



tion might affect the peak determination ( Comas- 



tri & Brusa 20081. The average black hole mass 



increases with increasing redshift (Merloni 2004), 
meaning that there are fewer low-mass black holes 
at high redshifts, and the average accretion rate 
decreases towards lower redshift. This observed 
anti-hierarchical behaviour of AGN evolution in- 
dicates that SMBHs and their host galaxies co- 
evolve. 

The underlying physics of the observed relation 
between the formation and evolution of galaxies 
and their central SMBHs is still unknown. Sev- 
eral efforts have been made to explain the ori- 
gin of this relation by using semi-analytical mod- 



Loeb 2003 ) and N-body SPH simulations ( Springel 
et al.||2005a|bl |Di Matteo et al.||2005| . Numeri 



cal simulations have shown that the gravitational 
tidal torques excited during major mergers lead 
to rapid inflows of gas into the centers of galax- 
ies ( Barnes fc Hernquist|1991 1996) which can be 
the mechanism to trigger quasar activity and star- 
bursts in galaxies. Observations of low-redshift 
quasars show a connection between galaxy mergers 
and quasar activity ( Heckman et al.|1984 Bahcall 
et al. 1996). In the local universe, ultraluminous 



infrared galaxies (ULIRGs) have bolometric lumi- 
nosities similar to bright quasars and are often in 



mergers (Sanders et al. 1986 Sanders fc Mirabel 



1996), and there is observational evidence that 



ULIRGs form the birthplaces of QSOs (Sanders 
eFaTlpSSl pan der Werf et al.|[20Tl| ). 



Di Matteo et al. (2005) argue that a merger 



generates a burst of star formation and leads to 
strong inflows that feed gas to the SMBH and 
thereby power the quasar. Hence, the energy re- 
leased by the quasar expels enough gas to inhibit 
further star formation and quenches the growth of 
the black hole. Also, Hopkins fc Quataert (2010) 



find that nuclear star formation is more tightly 
coupled to AGN activity than the global star for- 
mation rate of a galaxy. 

AGN, in which interstellar material is still 
rapidly accreting onto SMBHs, are key to the 
study of the Magorrian relation. The accretion of 
gas onto the central black hole yields a luminous 
source of X-ray, UV, and optical photons. Here, 
UV and X-ray radiation influence the chemistry 
of the accreting and star-forming gas and affect 
the thermodynamics of the ISM. Gas clouds in 
the inner regions of galaxies are exposed to radi- 
ation originating from the active nucleus, newly 
formed stars or both. The hot O and B stars 
that are formed in the starburst regions are ra- 
diating in the far-ultraviolet (FUV) wavelength 
range (6 < E < 13.6 eV), whereas hard X-rays 
(1 < E < 100 keV) are emanating from the nu- 
cleus when gas is falling in. When assuming that 
the energetics are dominated by either FUV or 
X-rays and effects by shocks can be neglected, the 
thermal balance and chemical structure are deter- 
mined by the radiation field and result in Photon- 



dominated Regions (PDRs, Tielens fc Hollenbach 
([1985])) or XDRs ([Maloney et al||1996[). 'AGN 



are thought to be partially obscured by dusty in- 
terstellar matter (ISM). X-rays (1-100 keV) have 
smaller absorption cross sections than UV pho- 
tons and, therefore, can penetrate large columns 
{N > 10^^ cm~^). Hence, they provide a probe 
of the innermost regions of AGNs. Furthermore, 
the main heating mechanism in X-ray dominated 
regions (XDRs) is photoionization (Coulomb heat- 
ing with thermal electrons), whereas in PDRs it 
is photoelectric heating. Hence, the heating efh- 
ciency in XDRs and PDRs is ~ 30% and 0.1-1%, 
respectively. We refer the interested reader to 



Meijerink & Spaans ( 2005 1 for a full introduction 
to the distinctions between PDRs and XDRs. 



Perez-Beaupuits et al. (2011 1 have post-processed 



the 3-D hydrodynamical models of an AGN torus 



(Wada et al. 2009) by using the chemical net- 
work of Meijerink & Spaans (20051. They esti- 
mated the X-ray flux emanating from the black 
hole and found that X-ray irradiation affects the 
thermodynamics of the AGN torus up to 100 
pc. They have calculated that the temperature 
of AGN gas, exposed to X-rays is a factor ^ 5 
higher than gas in a starburst of equal bolometric 
power. This has strong effects for star forming 
clouds, since the Jeans mass (Mj) scales with the 
temperature of the ambient gas as Mj c>c T^/^. 
Therefore, this might inhibit star formation in the 
central galactic regions or change the initial mass 
function [IMF, |Hocuk fc Spaans] ^OTo\ [20TT|)]. 



Moreover, Wada et al. (2009) studied the molec 



ular gas disks around AGN and have suggested 
that XDR physics may change the distribution of 
H2 around an AGN. Recent work by |Kim et al] 
(2011 ) focused on galaxy formation with feedback- 



regulated star formation and black hole growth. 
In their simulations they have taken into account 
the radiative feedback from the central MBH and 
found that radiative feedback from the MBH lo- 
cally suppresses star formation and self-regulates 
its growth. However, they use a monochromatic 
spectrum for X-ray photons and have a limited 
chemical network. 

In order to study the effects of X-ray irradi- 
ation from an accreting black hole in the cen- 
tral region of a galaxy, we perform 3-D cosmo- 
logical hydrodynamic simulations by implement- 



ing the XDR/PDR chemical network of Meijerink 
& Spaans (|2005) into Enzo. This work is the first 



X-ray irradiation by an AGN on ambient gas with 
a metallicity of zero and solar. We will study the 
Magorrian relation and derive observational diag- 
nostics to find the X-ray fingerprints in the high 
redshift universe in two subsequent papers. Here, 
as a proof of concept, we show that X-ray physics 
has been successfully included into Enzo. We fur- 
ther show that treating radiative feedback from 
the central MBH with a self-consistent XDR/PDR 
network and a polychromatic energy distribution 
is crucial to assess the relation between AGN and 
their host galaxies. In this study, we concentrate 
on the effects of X-rays on the multi-phase ISM 
near the central black hole. This will help us to 
better understand the conditions of star formation 
near an AGN and hence the interplay between stel- 
lar feedback and radiative feedback from the cen- 
tral BH. 

This paper is structured as follows. In Section 
2 we detail our implementation of the XDR code 
into Enzo and the treatment of the polychromatic 
spectrum. We further discuss the relevant cooling 
and heating processes in XDRs. In Section 3 we 
describe the set-up of the simulations. In Section 
4 we present the results and implications of our 
simulations. Finally, in Section 5 we discuss our 
results. We present some of the outputs of our 
chemical network in the Appendix. 

2. XDR IMPLEMENTATION 
2.1. Black Hole Radiation 

We set up our intrinsic specific X-ray flux (Fj) 
to have a spectral shape of the form 



^'(^)^^°(ilv 



cxp(— E/Ec) erg s ^ cm ^ 

(1) 

where Fg is the constant that determines the to- 
tal X-ray flux in the spectrum, E (> 1 keV) is 
the energy, a = 0.9 is the characteristic spectral 
index of the power law component of Seyfert 1 
galaxies ( Pounds et al.|1990, Madejski et al.|1995 



cV- 



part of a trilogy where we investigate the effects of 



Zdziarski et al.||1995[), and E^ = 10 keV is the 
cut-off energy ( Madejski et al.|[T995 ). In order to 
determine Eg, one needs an independent estimate 
of the bolometric X-ray flux Fx ■ We estimate the 
central Ex by assuming that only 10% of the total 
luminosity (Lf,o() is emitted in X-rays (Schleicher 
eFarilMol ), 



J>lkcV 

(2) 
where r is distance from the central black hole. We 
assume that the MBH has a radiative efhciency of 
e = 0.1 and thus has a luminosity L = e x Ledd = 
1.2 X 1037(e/0.1) M/Mq erg s'^. 

We set Mbh = 5 x 10^ M0, which gives a bolo- 
metric flux Fj^oi = 4 x 10^ erg s~^ cm~^ at our 
finest resolution of 3.6 pc. The UV part of the 
AGN spectrum is produced by the usual multi- 



color Shakura & Sunyaev (19731 accretion disk, 
with a viscosity parameter of a = 0.1. Here, a de- 
termines the efficiency of the transportation of an- 
gular momentum via a turbulent viscosity that is 
proportional to the local pressure in the disk. Our 
choice of a = 0.1 is consistent with observations 
( |Smak|1998[ |1999| |Dubus et al.|2001| and theoret- 



ical simulations of global accretion disks (Penna 
eraI][20T0l ). 



In our XDR models only X-rays with E > 1 keV 
are considered. Therefore, we take the integral of 
equation [l] over this energy range and equate it 
to the estimated bolometric flux to find Fg « 0.8 
erg s~^ cm~^ eV^^. Furthermore, with increas- 
ing distance from the AGN, the X-ray flux does 
not solely decrease by 1/r^ but also by attenua- 
tion along the line of sight. This attenuation is 
determined by the opacity, which is a function of 
energy and position 



T(E,r)=(7pa(E)NH(r), 



(3) 



where apa (E) is the photoelectric absorption cross 
section per hydrogen nucleus, and N//(r) is the to- 
tal column density of hydrogen between the cen- 
tral black hole and the position r. Taking this 
opacity effect into account we calculate the total 
X-ray flux for a given column density as 



'>ikcv IkeV 



)-ag-£/£.g-r(E,r)^^ (4) 



2.2. Radiative Feedback 

We use the Enzo radiation transport module 
Moray ( Wise fc Abel|20lT ) to calculate the trans- 
fer of X-ray radiation field produced by the black 



hole. Moray solves the radiative transfer equa- 
tion with ray tracing that is adaptive in spatial 
and angular coordinates. The ray normal direc- 
tions are determined by HEALPix (Hierarchical 
Equal Area isoLatitude Pixelation; |G6rski et al. 
2005). The rays originate from point sources and 



are progressively split when the angular sampling 
becomes too coarse, which can happen when the 
ray either propagates to larger radii or encounters 
a fine resolution AMR grid. 

The rays are traced through the grid in a typ- 
ical fashion (e.g., Abel et al. 1999) , in which we 
calculate the next cell boundary crossing. In each 
ray segment crossing a single cell, we solve the ra- 
diative transfer equation for a ray. 



IdP 
~c1H 



dP 

dr 



= -kP, 



(5) 



where P is the monochromatic photon flux associ- 
ated with the ray and k = «;(r, t) is the absorption 
coefficient. In the static case, its solution has a 
simple exponential analytic solution and the pho- 
ton flux of the ray is reduced by 



dP^Px (l-e^^). 



(6) 



where r is the optical depth (see above) across 
the ray segment. Because Moray equates the pho- 
toionization rate to the absorption rate, it is pho- 



ton conserving (Abel et al. 1999 Mellema et al 



2006). Thus, the photoionization and photoheat- 



ing rates are, respectively. 



"^ph — 



dP 



and 



"-absV^cllc'ip 



fcph (-Eph — Ei 



(7) 



(8) 



^ ph — "-ph ^-C'pli 

where V^cii is the cell volume, Pph is the photon 
energy and Ei is the ionization energy of the ab- 
sorbing material. In each cell, the resulting fcph 
and Fph values are the sum from all rays that cross 
that cell, which then are inputs into the chemistry 
rate equations and energy equation. 

2.3. Treatment of a General Input Spec- 
trum 



In the original release of Moray (Wise & Abel 



2011), radiation is discretized into energy bins. 



Since then, a more general approach of [Shapiro] 



et al.| ( |2004| and [Mellema et~aL] ([2006') has been 
implemented that can consider arbitrary (poly- 
chromatic) spectral shapes. In our work, we 
use a polychromatic energy distribution. In this 
method, the radiative transfer equation is numeri- 
cally solved before the simulation, giving a relative 
ionizing photon flux Ii, as a function of neutral hy- 
drogen column density Nh- The relative ionizing 
photon flux 



IANh) 



hv 



&v (9) 



for H I, He i, and He ii is computed and stored for 
300 column densities, equally log-spaced over the 



range Nh 



10 



12-25 



cm . fth is the ionization 



energy threshold for each absorbing species. Also, 
the average photon energy 



{E) ^ ■"'^'^ 



Qh,,UNH)du 

r/.(iVH)d. 



(10) 



is pre-calculated as a function of hydrogen column 
density Nh- The attenuation across each segment 
is thus dP = luiNH + dNn) - Iu{Nh), which is 
determined by interpolating from pre-calculated 
tables. These values of dP and {E) are used in 
Equations ([t]) and ([s]), resulting in a solution that 
retains the full spectral shape until the source ra- 
diation is fully attenuated. 

2.4. XDR physics 

The dominant parameter that drives the chem- 
ical and thermal structure of the gas in XDRs is 
the ratio of the X-ray energy deposition rate to gas 
density, Hjc/n. This occurs because the molecular 
destruction and heating rates per unit volume due 
to X-ray induced ionizations are proportional to 
nHjf, whereas the molecular formation and cool- 
ing rates are generally proportional to n^ times a 
rate coefficient ( Maloney et al.|[l996 ). 

The main heating mechanism in XDRs is pho- 
toionization. The rate at which photoionizations 
happen depends on the cross-section, which is a 
function of the frequency of the radiation and 
the properties of the atom. The absorption cross 
section of X-rays roughly scales with energy as 
^ 1/E^, which allows X-rays to penetrate deep 
into interstellar clouds. A photon with an energy 
of 1 keV penetrates a hydrogen column of about 
10^^ cm^^ for solar metallicity. 



UV radiation is absorbed by outer-shell elec- 
trons, whereas X-ray photons are absorbed by 
inner-shell electrons. When an electron from an 
inner shell (e.g., the K-shell) is ejected it leaves 
the atom behind in an excited state with an inner 
shell vacancy and more energy than is required to 
remove the least bound electron. The inner shell 
vacancy can be filled by an outer-shell electron 
(e.g., an L-shell electron). The excess energy then 
can be shed in one of two ways: by emitting pho- 
tons (fluorescence) or by ejecting outer electrons, 
which is referred to as the Auger effect. 

When inner-shell ionization is followed by 
Auger ionization, two high energy electrons are 
produced: the first from the primary ionization 
with an energy equal to the photon energy minus 
the binding energy of the ejected electron, and 
the second from Auger ionization which has an 
energy equal to the difference in binding energies 
of electrons in the inner and outer shell minus the 
binding energy of the outer shell electron. This 
can be hundreds or thousands of eV. In an ion- 
ized gas, these suprathermal electrons undergo 
frequent elastic collisions with thermal electrons 
and their kinetic energy is converted into heat. 

In a predominantly neutral gas, these suprather- 
mal electrons can coUisionally excite and/or ionize 
ions before interacting with free electrons. These 
secondary ionizations are more important for H, 
H2 and He than the primary ionizations, which is 
a consequence of the large primary photo-electron 
energies (Meijerink & Spaans 2005). The impor- 



tance of secondary ionization is determined by the 
ionization fraction x = n(H+)/n(H) since the prob- 
ability that a fast electron will share its energy 
with a free electron, compared to its probability 
of striking an atom or molecule, is proportional to 
this fraction. The ratio of secondary to primary 
ionizations is ^ 26, depending on the chemical 



composition of the gas (Dalgarno et al. 19991. If 



the ionization fraction of the gas is a; > 0.9, then 
most of the kinetic energy goes into heat through 
Coulomb interactions with the ambient thermal 
electrons. If a; < 0.9, then approximately 40% of 
the primary photo-electron energy goes into sec- 
ondary ionization and excitation (H I excitation, 
mainly Lya) whether the gas is atomic (Shull & 
van Steenberg||1985l iFurlanetto & Stoever||2010|) 



or molecular ( Voit||199T ). 

The heating efficiency of a molecular gas is 



larger than that of a mostly neutral atomic gas. 
When H2 is ionized by a fast electron and recom- 
bines dissociatively, about 10.9 eV of the ioniza- 
tion energy can go into kinetic energy. On the 
other hand, an Hj ion can react on with a hy- 
drogen molecule to produce Hg , which may then 
either recombine dissociatively or react with other 
species. Due to the strong exothermic nature of 
the re combination reaction, |Glassgold fc Langer| 
(1973) argue that for every Hj ion formed about 
8 eV goes into gas heating. 

In general, ro-vibrational H2 cooling is impor- 
tant in XDRs due to the high gas temperatures 
{Tg > 1000 K). The ro-vibrational bands of molec- 
ular hydrogen can be excited by the nonther- 
mal electrons produced through X-ray ionization 
and thermal collisions with e^, H and H2. Col- 
lisional de-excitation of vibrationally excited H2 
can be an important heating source when non- 
thermal electrons dominate the excitation (at low 
Hjjf/n). These excited H2 molecules can also en- 
hance chemical reactions with an activation bar- 
rier ( [Meijerink fc Spaans||2005 1 . 

XDRs are exposed to X-rays as well as FUV 
photons. Internally generated FUV photons 
are produced when energetic nonthermal elec- 
trons collide with atomic and molecular hydrogen. 
These collisions result in the emission of Lya and 
Lyman- Werner band photons, and significantly 
affect the chemistry of the X-ray irradiated gas 
clouds. If the electron fraction of the gas is x < 
0.01, then about 40% of the energy deposited by 
X-rays results in the production of these internal 
FUV photons. Thus, the dense interiors of X-ray 
irradiated clouds can chemically resemble pho- 
todissociation regions with significant abundances 
of neutral oxygen as well as singly ionized iron, 
silicon and carbon ( Maloney et al.|[T996 ). 

Most of the FUV photons that are locally pro- 
duced by the nonthermal electrons will be ab- 
sorbed by dust grains. The resulting tempera- 
ture of the dust grains is proportional to the lo- 
cally absorbed X-ray energy per hydrogen atom 
(Hx). Therefore, if the dust abundance is high, 
then there is less energy per dust particle and the 
average dust temperature T^ drops and is given 

by 



where Xd — 1.9 x 10~^ is the grain abundance 
(for solar metallicity) and Hx is in erg s~^ (Yan 
1997). The minimum grain size is set to B,jnin = 
lOA. In XDRs, the dust temperature for the same 
impinging flux by energy is lower than in PDRs 
(iMeijerink & Spaansl 120051). Also, in XDRs the 



gas heating efficiency is 10-50% whereas it is 0.1- 
1% in PDRs. Furthermore, the chemistry in XDRs 
is less stratified, with C+, C, CO co-existing over 
large columns, than in PDRs. 

At high temperatures (T > 5000 K) the gas 
cooling is dominated by coUisional excitation of 
Lya, and forbidden transitions of [O I] (AA 6300, 
6363 Aim), [C I] (AA 9823, 9850 ^m), [Fe II] (AA 
1.26, 1.64 ^m), and [Si II] (AA 6716, 6731^im). 

The cooling below a few 1000 K is dominated 
by the fine-structure lines of [01] 63 /im, [Sill] 35 
Atm, [CII] 158 /im, [CI] 269 and 609 /xm. Rota- 
tional and vibrational transitions of H2, H2O, and 
CO may also be important when these molecules 
are abundant. Cooling lines may become optically 
thick and such column density dependent transfer 
effects are included using a non-local escape prob- 
ability method ( jPoelman fc: Spaans 2005 ) . At high 
densities, and if the grain temperature is less than 
the gas temperature, gas-grain collisions can be 
a source of gas cooling and grain heating. For 
further details, we refer the interested reader to 



Meijerink fc Spaans (2005) 



2.5. Modifications to tlie 2005 XDR code 



The XDR code as described in [Meijerink fc| 
Spaans (2005) has been modified by including all 



the heavy elements up to iron with abundances 
> 10^^. As a result, the elements treated in the 
chemistry of the code were extended with Ne, Na, 
Ar and K. In the original XDR code, the doubly 
ionized states of C, N, O, S, and Fe were included. 
Here, we treat the singly and doubly ionized states 
of all elements, also He^+ . Mo st of the additional 
reactions were adopted from Adamkovics et al. 



Td = 1.5 X lO^Hx/xdf' K, 



(11) 



(2011), who give an elaborate description on the 
X-ray ionization of heavy elements. We followed 
their method to determine the secondary ioniza- 
tion rates of He^"*" , which were scaled to those for 
hydrogen, by using the peak electron impact cross 
sections as obtained from Tarawa & Kato (1987). 
The primary ionization rates were obtained from 



Verner fc Yakovlev ( 1995 1, where we assumed that 



no Auger ionization is possible, implying that only 



one electron is released per absorbed X-ray pho- 
ton. 

3. Simulations set-up 

In this work, we use the cosmological adaptive 
mesh refinement (AMR) code Enzo-2.oF] (Bryan & 
Norman 1997; O'Shea et al. 2004) that is modified 
to include XDR physics. We perform simulations 
in a three-dimensional periodic box with a side 
length of 3 h^^ Mpc, initialized aX z — 99. The 
size of the root grid is 128'^ with three nested sub- 
grids, each refined by a factor of two. The finest 
grid has an effective resolution of 1024^^ with a side 
length of 375 h~^ kpc. Refinement is restricted 
to the finest grid and occurs during the simula- 
tions whenever the baryonic matter, or dark mat- 
ter density, is greater than the mean density by 
a factor of four. The maximum level of refine- 
ment that is reached in the finest grid is 10, al- 
lowing us to have a resolution of 3.6 physical pc. 
We use Wilkinson Microwave Anisotropy Probe 
seven-year cosmological parameters (Komatsu et 
al. 2009), which have the following values: J^a = 
0.734, n^ = 0.266, n,, = 0.0449, erg = 0.81, and 
h = 0.701. Here, J^a is the vacuum energy, $!„ 
is the matter density, f2b is the baryon density, erg 
is the variance of random mass fiuctuations in a 
sphere of radius 8 ft.~^ Mpc, and h is the Hubble 
parameter in units of 100 km s^^ Mpc~^. 

In order to see the impact of X-ray radiation 
from the central MBH on the surrounding gas, we 
have performed two simulations. In these simula- 
tions we use the primordial chemical network of 
Enzo until z = 15, where the most massive halo 
with a mass of 2 x 10® M© forms. At z=15, we 
insert a MBH with a mass of 5 x 10^ M©, in or- 
der to be roughly consistent with the Magorrian 
relation, and consider its radiative feedback after- 
wards. In one of the runs, we start using the XDR 
tables compiled for solar metallicity gas (hereafter 
referred to XDRs), and in the other one, we keep 
using the primordial chemical network of Enzo 
(hereafter referred to XDRz). We compare a non- 
zero and zero metallicity run that both enjoy the 
full Enzo-XDR treatment. In both simulations, 
we employ the radiative transfer module, Moray, 
which uses a polychromatic energy distribution. 
Hence, the only difference between the two runs is 

^ http://enzo-project.org 



the metallicity of the ambient gas that is exposed 
to X-ray radiation from the accreting MBH. 

Our XDR chemical network consists of more 
than 170 species. We have constructed tables of 
XDR solutions for species abundances and tem- 
peratures for solar metallicity and for a wide range 
of X-ray flux Fx = lO^^^s - lO^^ erg cm'^ 
s~^, density n = 10-10^ cm^'^ and column den- 
sity Nh = 10^°-10^^ cm~^. This large parame- 
ter space enables us to model the ISM properties 
close to an AGN properly. Here, the XDR model 
assumes that the gas is instantaneously in local 
thermal equilibrium (LTE). In XDRs, the domi- 
nant heating mechanism is Coulomb heating and 
is transferred through collisions between electrons 
and gas. The slowest process of these two is the 
collisions. The particle collision time, which is the 
time scale on which heating occurs as electrons, 
generated by photoionization, thermalize with the 
ambient medium through elastic collisions and is 
given by l/{andv) in seconds. Here, n is the num- 
ber density in cm^"^, a is the cross-section in cm 
~^, and dv is the velocity in cm3 s~^. This cor- 
responds to t < 3 yr for densities in excess of 10^ 
cm^'^. Furthermore, at solar metallicity the cool- 
ing time kT/nA is < 30 yr for densities in excess 
of 10^ cm~'^. The time step of our simulations dt 
^ 30 yrs (on the order of 600 yrs). Hence, the 
LTE assumption holds. 

We use Enzo's 9 species (H, H+, H^, He, He+, 
He^"*", H2, H2", and e~) non-equilibrium chemical 



network for zero metallicity case (Abel et al.|1997 



Anninos et al. 19971. In the XDRs case, we run 



both the XDR routine and the chemical network 
of Enzo and compare the computed temperatures 
to determine if the cell is X-ray dominated. We 
take the highest value of the two found temper- 
atures and continue to iterate for the next step. 
By taking the maximum we basically divide the 
simulation box into XDR and non-XDR zones. Of 
course, we might overestimate the temperature in 
the grid when X-ray heating and non-X-ray heat- 
ing are comparable, but this pertains to a very 
small part of the grid given the deep penetration 
of X-rays into dense gas, where shocks are of mod- 
est importance. 

We choose solar metallicity for our XDR ta- 
bles since this is a good approximation for regions 
around SMBHs, even inside z ~ 6 AGNs (Jiang 
et al.|[2007). Moreover, attenuation of X-ray pho- 



tons in a non-zero metallicity ISM has never been 
done before for a high redshift accreting MBH, 
and here we show that we can treat the impact of 
X-rays on such non-zero metalhcity ambient gas. 
By choosing solar metaUicity, we seek to highhght 
the fuU range of effects that metals introduce in 
X-ray exposed gas. 

Furthermore, in our tables, we also store the 
abundances of all species that Enzo's chemical net- 
work is using, attached to the corresponding tem- 
perature values. For a given unattenuated flux, 
we create a table of attenuated fluxes as a func- 
tion of column density. Combined we thus have 
an XDR grid of models in n, Fx, and Njj (varying 
metallicity will be presented in paper II) that uses 
Moray to compute the full (chemical, thermal and 
hydrodynamic) response of X-ray exposed gas at 
non-zero metallicity. In the code, we calculate the 
column densities between the cell and MBH and 
incident flux with Moray, take the corresponding 
flux values from the pre-calculated tables and ap- 
ply the radial dependence. By using the flux, den- 
sity, and column density computed in Enzo, we 
find the corresponding temperature values and the 
species abundances from the pre-computed XDR 
grid tables for each cell and feed these back into 
Enzo. We stop the simulations after ^ 42 Myr at 
z=13.54 This is long enough for ambient gas in a 
radius of a few hundred pc to be affected by X- 
rays emanating from the 5 x 10'' Mq MBH and to 
find (if at all) an equilibrium again. 

In this work, we did not take into account the 
effects of momentum transfer from the ionizing ra- 
diation field. Addition of radiation pressure to cos- 
mological hydrodynamic simulations is not very 
well studied and hence we did not want to intro- 
duce another unknown parameter to our simula- 
tions. However, recently Wise et al. (20121 have 



performed the first cosmological simulations with 
radiation pressure that is calculated by solving the 
RT equation, and have shown that the main mech- 
anism for blowing away the gas from the central 
region of a galaxy is SNe. But one thing to keep 
in mind is that, in their simulations they reached 
maximum densities of ^ 100 cm^"^ and did not 



take into account dust. In his work, Draine (2011 ) 



we only concentrate on the X-ray effects and leave 
the investigation of the interplay between radia- 
tion pressure and X-ray irradiation for a follow up 
paper. 

For the analysis of our cosmological simulations 
we use yt, a cross-platform analysis toolkit written 
in Python ( [Turk et al.||2011 |. 



has made a systematic study on the effects of ra- 
diation pressure in a static, dusty H ii region and 
has shown that radiation pressure becomes impor- 
tant at high densities (n> lO'^ cm~^.) Hence, here 



4. Results & Implications 

In the simulations that we perform, as a proof of 
concept, we see significant differences between the 
solar metallicity and zero metallicity cases. Below 
we explain the physical processes that play a role 
in causing the differences in the ISM properties of 
the modelled halos. 

Shortly after we insert a MBH, we already note 
considerable differences between the two runs in 
the central region around the MBH. In Figure 
IT] we plot the density-temperature histograms, 
within a sphere of 500 pc diameter, for the XDRs 
(top row) and XDRz (bottom row) runs at red- 
shifts z = 14.95, 14.78, 14.54 and 13.54 from left 
to right. At z = 14.78, only 5.8 Myr after we turn 
on the central MBH, phase diagrams show that 
there is less gas present in the XDRz case than in 
the XDRg case. This difference in the gas mass, 
becomes even more pronounced after 12 Myr we 
turn on the MBH (z = 14.54). To estimate the 
amount of gas that is missing in the zero metal- 
licity case we compare the enclosed gas mass of 
both simulations as a function of radius at redshift 
z = 14.78 (left) and z = 14.54 (right) in Figure [5] 
We calculate the difference in enclosed gas mass 
(blue dot-dashed line) in the inner 200 pc to be 
8 X 10^ M0 between the XDRs and XDRz cases, 
after 5.8 Myr we turn on the MBH (z = 14.78). 

Furthermore, when we look at the 2D radial 
density and temperature profiles (see Figures Is] 
and [4]) of the XDRg and XDRz cases we see sig- 
nificant differences in the inner 100 pc while at 
larger radius (R > 500 pc) the differences become 
less pronounced, (see section 4.2 1. 

The missing mass in the inner 200 pc of the 
XDRz run can be explained as follows. In the 
XDRs case because of the metals. X-rays produced 
by the central MBH see a high opacity, about a 
factor of 3-30 higher than for zero metallicity gas. 
This is because metals like C, N, O, Si, and Fe 
have large cross sections, for inner shell absorp- 
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Fig. 3. — 2D profiles of density versus radius for the XDRs (top) and XDRz (bottom) at z = 14.95 (left), 
14.54 (middle), and 14.17(right). 



tions (see Figure 5.3 in Meijerink (20061 for E > 1 



keV and Morrison & McCammon (1983)), which 



lead to an "opacity wall" , if column densities ex- 
ceed 10^^ cm~^. Thus, the energy deposition rate 
into the medium close to MBH is high and the 
X-ray power is dissipated only locally. On top of 
that, the cooling ability of solar metallicity gas 
is high. Therefore, the cooling time is short be- 
yond the opacity wall. Combined, large column 
densities of 10^* cm~^ are reached in the central 
20 pc as metal enriched gas falls in and builds up 
sufficient column to shield itself from irradiation 
(see Figure l5|. However, in the XDRz case, in the 
absence of metals, the maximum column density 
of the bulk of the gas in the inner 20 pc (~ 10^° 
cm~^) is 3 orders of magnitude less than in the 
XDRs case aX i — 14.95 and the X-rays penetrate 
to larger columns. Thus, the energy deposition 
rate into the medium close to the MBH is signif- 



icantly less and X-ray power is dissipated glob- 
ally. This difference in column densities is shown 
in Figure [5j where we plot the 2D column density- 
radius profiles for both cases at redshift z = 14.95, 
14.54 and 14.17. In this plot it is easily seen that 
the gas column, viewed from the MBH, builds up 
much faster in the XDRs case due to the efficient 
accumulation of solar metallicity gas. Here, in the 
XDRs case there are no column densities shown at 
radius larger than 30 and 110 pc at z= 14.95 and 
14.54, respectively. This is due to the fact that we 
do not track the column density once the ray is 
terminated after it is almost completely absorbed. 
Therefore, in the XDRz case, the impinging X- 
ray flux onto the relatively low density gas that sits 
at larger radii is much higher. Hence, this forms 
an X — ray induced Hii region in the XDRz case 
already after only 5.8 Myr we turn on the central 
MBH. 
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Fig. 4. — 2D profiles of temperature versus radius for the XDRs (top) and XDRz (bottom) at z = 14.95 
(left), 14.54 (middle) and 14.17 (right). 



The H II ionization front will sweep through the 
surrounding medium until the recombination rate 
inside the H ii region balances the energy output 
of the MBH (see Figure [6] for the HI ionization 
fraction). As a consequence, in an attempt to es- 
tablish pressure balance, the newly formed H ii 
region will expand, thus driving large quantities 
of gas to larger radii. This also further lowers the 
column densities at small radii in the XDRz case. 
A back of the envelope calculation, where we com- 
pare the gas kinetic energy to the energy in X-rays, 
shows that the energy input from the black hole 
(Lfidd ^ 10'*^ erg s~^) is sufficient to drive ^ 10^ 
Mq gas out with an average speed of ^ 100 km 
s~^ in 6.5 Myr. 



4.1. Dynamics of the H II regions 

In Figure [7j we show the velocity magnitude 
at z = 14.78 (top) and 14.54 (bottom) for the 
XDRs (left column) and XDRz (right column) 
runs overlaid with velocity vectors showing the in- 
flow/outflow. Indeed, at z = 14.78 (after only 5.8 
Myr we turn on the MBH) we see an outflow in 
the XDRz case, at a speed of ~ 100 km s"'^, which 
is caused by the X-ray induced H II region. How- 
ever, in the XDRs case, due to the opacity wall, 
the X-rays are strongly attenuated. This leaves 
the gas at larger radii virtually unaffected. Due to 
the lack of a large scale H ii region in the XDRs 
run, the bulk of the gas keeps falling towards the 
MBH, increasing the densities around the black 
hole. One thing to keep in mind here is that our 
smallest cell size is 3.6 pc, which means that we 
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might be overestimating the metals' opacity effect 
and hence the energy deposition rate. It is pos- 
sible that higher resolution simulations will show 
the formation of an H ii region in the XDRs case 
as well, although at a much smaller spatial scale 
than in the XDRz run. In fact, at z = 14.17 (23 
Myr after we turn on the MBH), we see that X- 
rays start to penetrate much further (> 200 pc) in 
the XDRg case as well and drive an H ii region, but 
with a delay of 17 Myr with respect to XDRz case, 
which pushes apart the central opacity wall with 
a speed of ~ 200 km s~^. This is shown in Figure 
Is] where we plot the X-ray flux (top) and velocity 
magnitude (middle) slices along the x-axis, and 2D 
profiles of temperature vs radius for the XDRg run 
at z=14.17 (left column) and 13.54 (right column). 
The H II region forms after 23 Myr we turn on the 
MBH (z = 14.17), expands to a few kpc in 20 Myr 
at z = 13.54 (middle row) and drive the gas away 
from the central MBH (bottom row). 



To understand this delay, two effects are rele- 
vant. First, in the XDRz case, the initial H ii re- 
gion radius is larger because of the lower X-ray op- 
tical depth through zero metallicity gas. Second, 
for a local ionizing photon rate S^, the Stromgren 
radius scales as Rs oc {Si/n^Y^^, the recombina- 
tion time scales as t^ ex 1/n, and the ionization 
front velocity scales as V/ oc Si/nR'^. Because 
X-rays are absorbed in denser gas and at smaller 
radii in the XDRg run, the H ii region grows more 
slowly when compared to the XDRz case. 

4.2. Thermodynamics 

At small radius (R < 100 pc), as a consequence 
of high metal opacity, the attenuated X-rays heat 
up the gas in the XDRg case as shown in Figure 
|4j Therefore, the temperature in the inner 20 pc 
in the XDRg case increases from 10^ K to almost 
^ 10^ K. As a consequence of this the densities 
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drop from 10^^^ g cm^'^ to 10^^"* g cm^^ as shown 
in Figure [3) However, in the XDRz case, the den- 
sities at z = 14.95 in the inner 20 pc are 5 orders 
of magnitude less than in the XDRs case. The X- 
ray induced H 11 region sweeps away most of the 
gas from the inner 100 pc, keeping the densities 
low (10^^^ g cm~^) and temperatures around 10^ 
K at larger radii. 

At large radius (R > 500 pc), we do not see 
much difference in the multiphase ISM of the 
XDRg and XDRz runs. This is shown in Fig- 
ure|9]where we plot the density-temperature phase 
diagrams within a sphere of 1 kpc radius for the 
XDRg (top) and XDRz (bottom) runs at redshifts 
z = 14.95, 14.78, 14.54 and 13.54, from left to 
right. In a gas with modest ionization degree and 
weak X-ray radiation field the H^ route will drive 
the formation of H2 and HD, due to the avail- 
ability of free electrons. The presence of dust 



will boost the formation of H2 ( Cazaux & Spaans 



2009). However, an increase in H2 abundance will 



also increase the H2 ionization rate, and hence the 
gas heating rate, when X-rays are present. In our 
simulations, most of the gas has temperatures of 
a few thousand K, and at those thermally unsta- 
ble temperatures H2 is the dominant cooling and 
heating channel. Therefore, we do not see a funda- 
mental difference in the density and temperature 
profile of the ambient gas between the two runs. 

In the XDRs I'un we see that at redshifts z — 
14.78, as well as slightly later z — 14.54, high den- 
sity gas (> 10~^° g cm~'^) is more abundantly 
present at temperatures of > 500 K (see Figure 
[9]). There are two effects that play a role here. 
First of all, in the XDRg run the opacity is large 
enough to absorb any X-ray photon below 5-10 
keV. Also, we reach column densities of 10^"* cm^^ 
at small distances from the MBH, as shown in Fig- 



13 




5-10^ 1-10^ 1.5-10^ 

Velocity Magnitude (cm s"-*-) 



2-10^ 
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and XDRz (right column) runs at z = 14.78 (top) and 14.54 (bottom). 



ure[5] This explains the lower HI ionization rates 
in the XDRs case as shown in Figure [6J Along 
these large columns, the X-ray flux is almost com- 
pletely absorbed by the gas and thus the gas at 
larger distances is not much affected. Secondly, 
as mentioned before, solar metallicity gas can cool 
down and evolve to higher densities faster than 
zero metallicity gas, which was also found in the 



PDR studies of Aykutalp & Spaans (2011). Dur- 



ing the fast collapse, adiabatic heating raises the 
temperature of the gas in the XDRs case as oc r?!"^ 
while for X-ray heating this goes as oc n. However, 



in the XDRz run the maximum column density, 
within the central 200 pc, is only 10^" cm~^, and 
the X-rays are not fully attenuated in this cen- 
tral region. Thus, the incident X-ray flux at large 
distances is high, leading to a high heating rate 
and H2 formation rate. Furthermore, due to the 
long cooling times of the zero metallicity gas in the 
XDRz case, it takes longer to reach high densities 
and adiabatic heating is not significant. 
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5. Discussion 

As a proof of concept, we investigated the im- 
portance of metals in the vicinity of an Eddington- 
limited X-ray emitting MBH. The high tempera- 
tures (T > 10^ K, Figure ^ that are found in 
the inner regions around the MBH in the XDRs 
case might have severe consequences. These high 
temperatures can quench the gas accretion onto 
the MBH. Consequently, the accretion power de- 
creases and thus negative radiative feedback be- 
comes self-limiting. Furthermore, high opacity 
due to metals might aid star formation in some 
regions by shielding the high density gas against 
X-ray irradiation and allow gas to cool and col- 
lapse. On the other hand, in metal-poor regions, 
if an X-ray induced H ii region forms, the lat- 
ter can either quench star formation by removing 
gas or induce star formation by gas compression: 
As shown in our zero metallicity simulation, the 
X-ray driven H ii region sweeps out a significant 
amount of mass which means that the growth of 



the black hole will be reduced. Moreover, this 
would also quench the star formation in the in- 
ner regions around the black hole. But, the H ii 
ionization front might also induce star formation 
by compressing gas where otherwise stars would 
not form in gravitationally stable gas. Also, the 
redistribution of gas will cause a much flatter den- 
sity profile in the central region than in the high 
metallicity case. 

In this work, we showed that the difference be- 
tween zero and solar metallicity gas irradiated by 
X-rays is very large and that in order to study the 
Magorrian relation a self-consistent treatment of 
the effects of metals on X-ray physics should be 
included. Here, we only concentrate on the X-ray 
effects and did not take into account the radiation 
pressure. As discussed earlier, the effects of mo- 
mentum transfer from the ionising radiation field 
might further help to sweep away the gas from 
the central region of the MBH when densities of 
n> 10^ cm~^ are reached. Also, in our simulations 
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we compare the computed temperatures from the 
XDR routine and the chemical network of Enzo 
and take the highest value of the two found tem- 
peratures. By taking the maximum we basically 
divide the simulation box into XDR and non-XDR 
zones. Of course, we might overestimate the tem- 
perature in the grid when X-ray heating and non- 
X-ray heating are comparable, but this pertains to 
a very small part of the grid given the deep pen- 
etration of X-rays into dense gas, where shocks 
are of modest importance. In a follow-up paper, 
we will include a star formation recipe, a Bondi 
Hoyle accretion rate for the growth of the MBH 
and H2 self-shielding. 

Appendix 

In the following, wc show a few models from the 
grid of XDR models. These are merely intended to 
illustrate the dynamic range and variety in XDR 
physics that Enzo has been augmented with. A 
copy of the electronic tables will be made available 
upon request to the authors. 
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Fig. 10.— Abundances of e", H, H+, H", H2, H+, He, He+, Hc2+, C, C+, 0^+, O, 0+, 0^+, CO and 
temperature for column densities of 10^'^ cm^^ (left column) and 3.16 x 10^^ cm^^ (right column) for a 
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Fig. 11.— Abundances of e", H, H+, H", H2, H+, He, Hc+, Hc2+, C, C+, C^+ , O, 0+, 0^+, CO and 
temperature for column densities of 10^° cm^^ (left column) and 3.16 x lO^'' cm^^ (right column) for a 
X-ray flux of log F = 2.75 - 3.5 erg s~^ cm"^. 
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